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Abstract 

We consider the use of preconditioning methods to accelerate the convergence to a steady 
state for both the incompressible and compressible fluid dynamic equations. We also consider 
the relation between them for both the continuous problem and the finite difference approxi- 
mation. The analysis relies on the inviscid equations. The preconditioning consists of a matrix 
multiplying the time derivatives. Hence, the steady state of the preconditioned system is the 
same as the steady state of the original system. For finite difference methods the preconditioning 
can change and improve the steady state solutions. An application to flow around an airfoil is 
presented. 
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1 Introduction 


Seymour Parter has considered preconditioning methods for numerical approximations to elliptic 
partial differential equations. As an extension of his ideas we shall consider similar techniques for the 
fluid dynamic equations. Much effort has been expended to solve the compressible steady state fluid 
equations for a large range of Mach numbers. A standard way of solving the steady state equations 
is to march the time dependent equations until a steady state is reached. Since the transient is 
not of any interest one can use acceleration techniques which might destroy the time accuracy but 
enables one to reach the steady state faster. Such methods can be considered as preconditionings 
to accelerate the convergence to a steady state. For the incompresible equations the continuity 
equation does not contain any time derivatives. To overcome this difficulty, Chorin [2] added an 
artificial time derivative of the pressure to the continuity equation together with a multiplicative 
variable, /? . With this artificial term the resultant scheme is a symmetric hyperbolic system for the 
inviscid terms. Thus, the system is well posed and and numerical method for hyperbolic systems 
can be used to advance this system in time.. The free parameter /3 is then chosen to reach the steady 
state quickly. Later Turkel ([15], [16]) extended this concept by adding a pressure time derivative 
also to the momentum equations. The resulting system after preconditioning is no longer symmetric 
but can be symmetrized by a change of variables. 

It is well known, that it is difficult to solve the compressible equations for low Mach numbers. 
For an explicit scheme this is easily seen by inspecting the time steps. For stability, the time step 
must be chosen inversely proportional to the largest eigenvalue of the system which, for slow flows, 
is approximately the speed of sound, c. However, other waves are convected at the fluid speed, u , 
which is much slower. Hence, these waves don’t change very much over a time step. Thus, thousands 
of time steps may be required to reach a steady state. Should one try a multigrid acceleration one 
finds that the same disparity in wave speeds slows down the multigrid acceleration. With an implicit 
method an ADI factorization is usually used so that one can easily invert the implicit factors. The 
use of ADI introduces factorization errors which again slows down the convergence rate when there 
are wave speeds of very different magnitudes [12] . 

For small Mach numbers it can be shown ([5], [7], [9]) that the incompressible equations ap- 
proximate the compressible equations. Hence, one needs to justify the computational use of the 
compressible equations for low Mach flows. We present several reasons why one would still use the 
compressible equations even though the Mach number of the flow is small. 

• There are many highly efficient compressible codes available that could be used for such 
problems especially in complicated geometries. 

• For low speed aerodynamic problems at high angle of attack most of the of the flow consists 
of a low Mach number flow. However, there are localized regions containing shocks. 

• In many problems thermal effects are important and the energy equation is coupled to the 
other equations. Then, the compressible equations must be used even for low Mach number 
flows. 

Therefore, one wants to change the transient nature of the system to remove this disparity of the 
wave speeds. Based on an analogy with conjugate gradient methods such methods were called [15] 
preconditioned methods since the object is to reduce the condition number of the matrix. Another 
approach, in one dimension, is to diagonalize the matrix of the inviscid term. One can then use a 
different time step for each equation, or wave. Upon returning to the original variables one finds 
that this is equivalent to multiplying the time derivatives by a matrix. Hence, this same approach 
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is named characteristic time stepping in [17]. In multidimensions one can no longer completely 
decouple the waves and so the characteristic time stepping is only an approximation. 

Thus, for both the incompressible and compressible equations we will consider systems of the 
form 


w t + }x + 9y — 0. 

This system is written in conservation form though for some applications this is not necessary. Our 
analysis will be based on the linearized equations so that the conservation form does not appear 
in the analysis though it does appear in the final numerical approximation. This system is now 
replaced by 

P _1 Wi + fx + 9y = o, 

or in linearized form 

P” 1 !^ + Aw x + Bw y = 0, (1) 

with A and B constant matrices. 

In order for this system to be equivalent to the original system, in the steady state, we demand 
that P -1 have an inverse. This only need be true in the flow regime under consideration. We shall 
see later that frequently P is singular at stagnation points and also along sonic lines. Thus, we 
will temporarily consider strictly subsonic flow without a stagnation point. For transonic flow it is 
necessary to smooth out the singularity in a neighborhood of the sonic line. We also assume that 
the Jacobian matrices A = and B = are simultaneously symmetrizable. In terms of the 
‘symmetrizing’ variables we also demand that P be positive definite. We shall show later, in detail, 
that it does not matter which set of dependent variables are used to develop the preconditioner. 
One can transform between any two sets of variables. Popular choices are two out of density, 
pressure, enthalpy, entropy or temperature in addition to the velocity components. Thus, when 
we are finished we wiE analyze a system which is similar to (1), where the matrices A and B are 
symmetric and P is both symmetric and positive definite. Such systems are known as symmetric 
hyperbolic systems. One can then multiply this system by w and integrate by parts to get estimates 
for the integral of wf 9 i.e. energy estimates. These estimates can then be used to show that the 
system is weU posed (see e.g. [5]). We stress that if P is not positive then we may change the 
physics of the problem. For example, if P = —I then we have reversed the time direction and must 
therefore change aE the boundary conditions. Hence, to be sure that the system is well posed with 
the original type of boundary conditions we shaE only consider the symmetric hyperboEc system. 

With these assumptions the steady state solutions of the two systems are the same. Assuming 
the steady state has a unique solution, it does not matter which system we march to a steady state. 
We shaE later see that for the finite difference approximations the steady state solutions are not 
necessarily the same and usuaEy the preconditioned system leads to a better behaved steady state. 

2 Incompressible equations 

We first consider the incompressible inviscid equations in primitive variables. 

U X + Vy - 0 

u t + uu x + vu y +p x = 0 

V t + UV X + VVy+Py - 0 
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We consider generalizations of Chorin’s pseudo-compressibility method [2]. Using the precondition- 
ing suggested in [15] (with a = 1) we have 


or in conservation form 


pPt +U X + Vy 


u 


—Pt + U t + UU X + VUy + p x 
JpPt + V t + UV X + l)Vy + Py 


0 

0 

0 


P 2 


Pt + Ux 4" v y — 0 


j^Pt + + (w 2 + p)x 4- (u«)v = 0 

-jpPt + V t + (uv) x + ( v 2 + p)y = 0 


We can also write (2) in matrix form using 

( IIP 2 0 0 


P 1 = 


v/P 2 0 1 



( P 2 

0 

°\ 

1 9 

P = -« 

1 

0 

/ 

V- 

0 
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( 2 ) 
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1 IP 2 

0 

0\ 


(p\ I 

f 0 
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0 

u/P 2 

1 

0 


u 4- 

1 

u 

0 

v/P 2 

0 

U 


\ v It ' 


0 

u 



Multiplying by P we rewrite this as 


We also define 



w t + P Aw x + PBwy = 0. 


D = u>\ A + u>iB -1<wi,u>2<1 


= 0 


(3) 


where u>\ , u-j, are the Fourier transform variables in the x and y directions respectively. The speeds 
of the waves are now governed by the roots of det(XI — P — P Bu<i) = 0 or equivalently 
det{ AP -1 — Au\ — Bu 2 ) = 0. Let 


q = uuii + vu> 2 - 


Then the eigenvalues of PD are 


do = q (4) 

d± = ±P 
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and so the ‘acoustic’ speed is isotropic. 

The spatial derivatives involve symmetric matrices, i.e. D is a symmetric matrix but P is not 
symmetric. Thus, while the original system was symmetric hyperbolic the preconditioned system 
is no longer symmetric. In [15] it is shown that as long as 

0 l > ( u 2 + v 2 ) 

then the equations can be symmetrized. On the other hand the eigenvalues are most equalized if 
/3 2 = (u 2 + v 2 ) [15]. Hence, we wish to choose /? 2 slightly larger than u 2 + v 2 . However, numerous 
calculations verify, that in general, a constant /? is the best for the convergence rate. The reasons 
for this are not clear. 

We wish to stress that /? has the dimensions of a speed. Therefore, /3 cannot be a universal 
constant. There are papers that claim that /3 = 1 or (3 = 2.5 are optimal. Such claims cannot be 
true in general. It is simple to see that if one nondimensionalizes the equation then fi gets divided 
by a reference velocity. Hence, the optimal ‘constant’ (3 depends on the dimensionalization of the 
problem and in particular depends on the inflow conditions. In many calculations the inflow mass 
flux is equal to 1 or else p + (w 2 4- v 2 )/2 = 1. Such conditions will give an optimal /? close to one. 
However, if one chooses the incoming mass flux as ten then the optimal (3 will be larger. 

We next define the Bernoulli function 

H = p+(u 2 + v 2 )/2. 

Bernoulli’s theorem states that for steady inviscid flow H is constant along streamlines. We now 
multiply the second equation of (2) by u and the third equation of (2) by v and add these two 
equations. If 0 1 = u 2 + v 2 , the result is 


H t + uH x + vH y = 0. (5) 

Thus, by altering the time dependence of the equations we have constructed a new equation in 
which H is convected along streamlines. Furthermore, if H is a uniform constant both initially and 
at inflow then H will remain constant for all time. On the numerical level this will usually not be 
true because of the introduction of an artificial viscosity or because of upwinding. For viscous flow, 
(5) is replaced by 

H t + uH x + vH y = ^-(uAu + t?Au) 

We note that these relationships for H follow from the momentum equations and do not depend 
on the form of the continuity equation. Hence, we consider the following generalization of (2) 


1 rT 

— Pi + aH t + u x + v y — 0 
an 

-pPt + u t + uu x + vu y + p x = 0 

av 

JpPt +V t + UV X + VVy+Py = 0 


where, a is a free parameter. The eigenvalues of PD are independent of the parameter a and are 
given by (4). For a = 0,a = 1 we recover our original scheme. For a = — 1 the time derivative of 
the pressure no longer appears in the continuity equation. For general f) we have 


4 



p 



(a + 1) au av 
au p 2 0 

av 0 ft 2 


) 



P 2 

- au 
- av 


-au 


1 + a-*f£ 

aotuv 


-av 

aauv 


1 + a — 


\otu* 

W 


\ 

/ 


where d = 1 + a — aa^-jp?- and we require that d > 0. If P 2 = u 2 + v 2 and a = 1 then 


( u 2 + v 2 -au —av \ 

- 1 + 3^ tffc. 

In [16] an analogy to the symmetric preconditioning of van Leer, Lee and Roe was constructed for 
the incompressible equations. If we choose a = l,a = 1 we get this preconditioning of van Leer 
et.al., i.e. P is symmetric . 

These examples show that the preconditioning is not unique. If fact, since the determinant of the 
transpose of a matrix is equal to the determinant of the original matrix it follows that the transpose 
of P is also a preconditioner with the same eigenvalues for the preconditioned system. In general, 
these various systems will have similar eigenvalues but different eigenvectors for the preconditioned 
system. Numerous calculations show that the system given by P in (2) is more robust and converges 
faster than that with the transpose preconditioner. This shows that it is not sufficient to consider 
just the eigenvalues but that the eigenvectors are also of importance. However, even when P is 
symmetric PD is not symmetric and so the eigenvectors of the preconditioned system do not form 
an orthogonal basis. 

We next examine some general form that the preconditioner can have. For this analysis it is 
easier to use streamwise coordinates as suggested in [17] and so v = 0. Let u m be some normalization 
of the velocity components, then 

( 0 u. 0 \ / 0 0 u* \ 

u, ti 0 , fl = I 0 0 0 j 

0 0 u / \ U. 0 0/ 


Then the "convective” eigenvector for the non-preconditioned system is 


0 

U>2 

-Wi 


The "acoustic” eigenvectors are given by 


-uu»i +'^/(uwi ) 2 +4uJ 
2 

U*U>1 

li*W2 


— VCJl --\/(txwi) 2 +4u? 
2 

U+U)\ 

U+ U>2 


5 


We now consider preconditioners of the form 


[a 6 0 \ 
P= c d 0 . 

\0 0 1 ) 


( 6 ) 


Let D = u\A + u 2 B , — 1. We want the eigenvalues of PD to be W\u, ±u. This gives 

us three relations for the four unknowns: 



( 6 + c)u* + du = 0 
u 2 d — bcu 2 = u 2 

The values suggested in [15] are b = 0,c = = 1 while the values suggested in [17] are 

b = c = — ^,d = 2 . We next present the eigenvectors of PD in terms of the elements of P. We 
exclude the case u 2 = 0, = ±1 as in this case PD has a double eigenvalue u and the eigensystem 

completely changes. Then the ’’convective” eigenvector is 



The "acoustic” eigenvectors are given by 


t i»(a + bu%) - u,(6 + c)u>i 
(u»6a>i + u)u>2 


£-buu>l+^-bvu , 
u.(a + fcwj) + u,(b + c)u>] 
(tt.ftwi — «)W2 


We note that the convective eigenvector is the same as before the preconditioning for the choice 
6 = 0. The two acoustic eigenvectors are orthogonal to each other if we choose 6 = 0 and c 2 = 
« 2 (i-£) 

— . This is similar, but not identical, to the choice suggested in [15]. There is no way to make 
the convective eigenvector normal to both acoustic eigenvectors for preconditioners of the form (6). 


3 Compressible equations 

The time dependent Euler equations can be written as 


Pt + up x + vp y + pa 2 {u x + t> y ) = 0 

Px 

u t + uu x + vUy + — = 0 (7) 

Pu 

v t + UV X + VVy + -2- = 0 

p 

S t + uS x + vS y = 0 

where a is the speed of sound given by a 2 — 

The form of this system is unchanged if we nondimensionalize the equations. From now on 
we shall assume that u,v,p,p are nondimensional quantities where the dimensional variables are 
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nondimensionalized by u.,p.,p. with p* = p.u 2 . Following [5] we define e = . If the fluid is 

isentropic then 

( 8 ) 


£_ 
P- 7e 2 


and 


_ £ 




( 9 ) 


Hence, as e goes to zero the speed of sound, a, goes to infinity and so the first equation in (7) 
reduces to u x + v v = 0. 

It was pointed out in ([15], [16]) that these equations can be symmetrized by using as the 

independent variable rather than dp . Hence, we define a new variable <f> by d<f> = For isentropic 
flow both p and a are functions only of the density and so using (8, 9) this can be integrated 

explicitly. This gives <f> = As the Mach number goes to zero <f> tends to infinity and therefore, 

Gustafsson and Stoor [5] subtract a constant and define 


<t> = 


__ p ¥-1 




This amounts to specifying the constant in the integration of d<f> from dp. They then prove, using 
energy methods, that for the linearized equations 

, dpincompretsible 

te 

Since p — ► 1 and using the definition of d4> this is equivalent to 

dp compressible * dpincompressible • 

We consider preconditionings that are a generalization of (3) 

/ £ 0 0 0 W * \ (u a 0 »\ (*\ 


( 10 ) 


^10 0 
^010 
Vo 001/ 


pa 

du 
dv 
\ dS ) 


a u 0 0 

0 0 u 0 

V 0 0 0 u ) 


du 
dv 

\dS ) 


O 

© 

o 


& \ 
pa 

0 v 0 0 


du 

a 0 v 0 


dv 

o 

o 

o 


\ dS J 


= 0 


The nonpreconditioned case corresponds to f3 2 — a 2 , o: = 0. Let 

q = uu\ + vu >2 

then the eigenvalues of PD are given by 

do = q ( double ) 

d± = 1/2 P 


(1 - a + /3 2 /a 2 )q ± y/((l - a + /3 2 /a 2 )V + 4(1 - g 2 /a 2 )/? 2 ] 

If we consider the special case a = 1 + /3 2 /a 2 we find that the ‘acoustic’ eigenvalues are given by 

d± = \/(l - q 2 /a 2 )0 2 


7 


Hence, these eigenvalues are isotropic in the limit of M going to zero. However, this eigenvalue 
vanishes at the sonic line and so the matrix is singular. In general, if we demand that the acoustic 
eigenvalues be isotropic then we have a singularity at the sonic line where the eigenvalues cannot 
be isotropic. The two ways out of this difficulty are either to smooth the formulas near the singular 
line or else to give up on isotropy. This difficulty is not a property of the preconditioning just 
presented but applies equally to all preconditioners. 

We now consider the system (7) in conservation form. 

Pt + (PU) X + ( pv)y = 0 

(pu) t + ( pu 2 + p) x + (puv) y = 0 

(pv) t + (puv) x + (pv 2 + p)y = 0 

E t + (pHu) x + ( pHv) y = 0 

where 

p u 2 + v 2 
E = - J — + — - — 

7-1 2 

E + p a 2 u 2 + v 2 

P ~ 1- 1 + 2 

Note that the Bernoulli function H is not identical with H for the incompressible equations. How- 
ever, we again have that for steady in viscid flow H is constant along stream lines. We now pre- 
condition the density and the energy equations in the following consistent manner. Let ^ be any 
variable we choose. Then we consider 

+ {py)x + ( pv)y = o 
{ipH)t + ( pHu) x + (pHv)y = 0 

Manipulating these equations gives 


H t + uH x + vH y = 0 


i.e. the total enthalpy, H, is simply convected in time along streamlines as we obtained for H in the 
incompressible case. It is interesting to observe that in the incompressible case we achieved this 
by preconditioning only the momentum equations while for the compressible flow we achieve this 
by preconditioning the continuity and energy equation. Of course, for isentropic flow the energy 
equation is not independent of the other equations and the result is not surprising. 

For the finite difference equation this implies that the artificial viscosity for the continuity 
equation should be based on ^ smd for the energy equation on ipH . If we choose = /?, i.e. no 
preconditioning for the continuity equation then we have the same artificial viscosity as suggested 
in [6] but with a different variable being advanced in time. If we choose ^ = p then both the 
continuity and energy equations are preconditioned. 

We next present the van Leer-Lee- Roe preconditioning for general non-aligned flow in (^, du, dv, dS) 
variables [17]. 


Pv = 


V 


— ■jpu/a 

— fiv/a 

0 


— pju/a 




(i? + l)^r 


uv 

+V 2 


— ~@i v l a o 

+ 0 

(^ + 0 

0 1 
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At the sonic line /? = 0 and r = 0 and the preconditioning matrix becomes singular. This 
preconditioning is not unique even if one only considers symmetric preconditioners. In both these 
examples the preconditioner was constructed based on using (p,u,v,S) as the dependent variables. 
The reason for this choice is that the matrices are symmetric which this choice. However, if 
another choice of variables is more appropriate tht introduces no difficulties. Thus, for example 
[1] recommends the use of (p, u, v,T) variables for the Navier-Stokes equations. Given two sets of 
dependent variables w and W let W w be the Jacobian matrix Then, we have dW — W w dw. So 
we can go between any sets of primitive variables or between primitive variables and conservation 
variables. In particular since the equations are solved in conservation variables we have several 
ways of going from the primitive variable preconditioner to a conservation variable preconditioner. 
Thus, the choice of variables used in constructing the preconditioner is dictated by mathematical 
or physical reasoning and then the preconditioner can be transformed to any other set of variables. 

• Construct the preconditioner matrix for the conservation variables. If W are the conservative 
variables and w the primitive variables then P conservative — (W^,) ^ P primitive ( Wyj ) ■ Details 
of the matrix Jacobians between various sets of variables are given in the appendix. 

• We calculate the residual dW in conservative variables. We then transform dW to dw as 
before. Next we multiply by P and finally transform back to conservative variables dW and 
update the solution. This is algebraically equivalent to the first option but requires three 
matrix multiplies instead of one. However, it offers more flexibility. 

• Similar to the previous suggestion we calculate the residual dW and transform to conservative 
variables dw and the multiply by P. At this stage we update the primitive variables w. We 
then use the nonlinear relations to construct W from w. This approach has advantages if the 
boundary conditions are given in terms of the primitive variables (p or T) and so they can 
be specified exactly and not approximately. 

If the residual dW is kept from the conservation form but the time derivative Wt is replaced 
by the time derivative of other variables, Wt this is linearly equivalent to preconditioning by 
the matrix P -1 = §|£. 

These methods are all equivalent for linear systems and the difference between them is mainly 
one of convenience. However, we shall next see that for the difference approximation these ap- 
proaches are not equivalent. 

4 Difference Equations 

Until now the entire analysis has been based on the partial differential equation. We now make 
some remarks on important points for any numerical approximation of this system. 

• For an upwind difference scheme based on a Riemann solver this Riemann solver should be 
for the preconditioned system and not the original scheme. In [3] plots are shown to illustrate 
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the greatly improved accuracy for low Mach number flows when the Riemann solver is based 
on the preconditioning. 

• For central difference schemes there is a need to add an artificial viscosity. Accuracy is 
improved for low Mach number flows if the preconditioner is applied only to the physical 
convective and viscous terms but xiot to the artificial viscosity. Volpe [19] shows that the 
accuracy of the original system deteriorates as the Mach number is reduced. The use of a 
matrix artificial dissipation ([14]) should be based on the preconditioned equations as in the 
upwind difference scheme. Upwind schemes without preconditioning tend to have difficulties 
with accuracy for low Mach flows [3]. 

Hence, both for upwind and central difference schemes the Riemann solver or artificial viscos- 
ity should be based on P _1 |PA| and not \A\. i.e. in one dimension solve Wf+P f x — ( \PA\w x ) x 
. For a scalar artificial viscosity |PA| is replaced by the spectral radius of PA or equivalently 
the time step associated with the preconditioned matrix. This is equivalent to not multiplying 
the artificial viscosity by P. 

• For a central difference scheme with a scalar artificial viscosity the artificial viscosity is of 
the form of a high order difference of the same quantity as is advanced in time. Thus, the 
continuity equation is solved for the density and so the artificial viscosity is a difference of the 
density. Similarly, for the momentum equations. For, the energy equation one can base the 
artificial viscosity on the energy. Alternatively it can be based on the total enthalpy which 
guarantees, for inviscid flow, that the total enthalpy is constant in the steady state [6]. When 
preconditioning the system one of the alternatives described above was to replace the time 
derivative of the conservative quantities with the time derivative of other variables. This, 
implies that the artificial viscosity should also be changed. Thus, if the continuity equation 
is updated for the pressure rather than the density, then the artificial viscosity should be 
based on the pressure. This is physically more reasonable for low speed flow since the density 
is almost constant and so will not contribute any reasonable viscosity. Furthermore, using 
a viscosity _i_n the continuity equation based on the pressure mimics what was done for the 
incompressible equations. This allows the low speed compressible equations to replicate the 
results of the incompressible equations on the finite difference level. This will be discussed in 
more detail in the following sections. 

• When using characteristics in the boundary conditions these should be based on the charac- 
teristics of the modified system and not the physical system. 

• When using multigrid it is better to transfer the residuals based on the preconditioned system 
to the next grid since these residuals are more balanced than the physical residuals. 

Preconditioning is even more important when using multigrid than with an explicit scheme. 
With the original system the disparity of the eigenvalues greatly affects the smoothing rates 
of the slow components and so slows down the multigrid method, [10]. 

• As indicated above there are accuracy diffi culies at low Mach numbers [19]. Some of these 

can be alleviated by preconditioning the dissipation terms. For very small Mach numbers 
there is also a difficulty with roundoff errors as — ► oo. Several people have suggested 

subtracting out a constant pressure from the dynamic pressure. A more detailed analysis 
[4] suggests replacing the pressure p by p where p = and e is a representative Mach 

number. 
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We conclude from the above remarks that the steady state solution of the preconditioned sys- 
tem may be different from that of the physical system. Thus, on the finite difference level the 
preconditioning can improve the accuracy as well as the convergence rate. 

In the previous section we stated that it is not important if one updates a different set of vari- 
ables or else uses the conservation variables and compensates with preconditioning by a matrix 
multiplication. However, numerically for very small Mach numbers the entries in the precondition- 
ing matrix can become very large or small. Hence, it can be advantageous to update the pressure 
or temperature directly rather than using a matrix multiply for preconditioning. 

5 Convergence 

We have previously quoted several papers ([5], [7]. [9]), that prove the convergence of the compress- 
ible equations to the incompressible equations, for isentropic flow, as the Mach number goes to 
zero. For nonisentropic flow there are no formal proofs. However, it is clear that for viscous flows 
that the boundary condition on the temperature, adiabatic or isothermal is very important, see 

[ii]- 

All these results refer to the time dependent physical equations. Once preconditioning is intro- 
duced time accuracy is lost and one can only discuss convergence of the steady state solutions. In 
this case one would hope that the time dependent preconditioned compressible equations converge 
to some time dependent preconditioning of the incompressible equations. In addition, one would 
also want this to be true on the numerical level. Thus, one would want to solve the preconditioned 
compressible equations by some numerical technique, on a fixed mesh and compare that with the 
solution of the incompressible equations on the same mesh. Mathematically, we have two limit pro- 
cesses occuring: the Mach number going to zero and the mesh size going to zero. These two limits 
need not commute. If one first converges the mesh size and then the Mach number it is equivalent 
to the convergence proofs for the analytic case. The more interesting problem is to converge the 
Mach number and then converge the mesh, i.e. we use a fix mesh as the Mach number is reduced. 

In particular this requires a careful study of the viscosities introduced by the scheme. We first 
consider an upwinding scheme. For the compressible case we have already noted that the Riemann 
solver should depend on the preconditioned problem. One would then need to show that this Rie- 
mann problem converges to a Riemann problem for some preconditioning of the incompressible 
equations. We next consider a central difference scheme with a scalar viscosity. In this case a high 
order even difference of some quantity is added separately to each equation, e.g. for the incom- 
pressible equations: pressure for the continuity equation, u and v for the momentum equations. For 
the compressible equations one normally adds a density difference to the continuity equation. In 
such a case it is obvious that the numerical scheme for the compressible equations cannot converge 
to the numerical scheme for the incompresssible equations. Furthermore, for low Mach number 
flows the density is almost constant and so the higher order difference of the density does not add 
much of a viscosity to the continuity equation. As such, we conclude that the artificial viscosity 
for the compressible continuity equation should be based on pressure and not density (at least for 
low Mach numbers). 

We shall examine the convergence a little more closely. By convergence of the compressible 
equations to the incompressible equations we are merely verifying what happens to the difference 
equations as the Mach number goes to zero. The convergence of the solution of the numerical 
approximation to the preconditionined compressible equations to the numerical solution of the 
incompressible equation is more difficult. However, we shall see that for the numerical solution the 
convergence of the difference equation is nontrivial and depends on the preconditioning. For this 
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purpose we shall only consider a central difference approximation together with a scalar artificial 
viscosity for the nondimensionalized preconditioned inviscid equations. 

For the incompressible equations in nonconservative form we consider the preconditioned system 


Pt + "h tfy) — h [(^lPut)i + (KiPyyy^y] 

■j^Pt + U t + UU X + VUy + p x = h, 3 [(K\U XXX ) X + (K 2 Uyyy)y] ( 11 ) 

Pt + V t + UV X + VVy +Py = k 3 [(KlV xxx ) x + (K 2 Vyyy)y] 

where each space derivative is approximated by a central difference with spacing h in each direction. 
The time derivatives are replaced by a multi-stage scheme. K\,K 2 are the largest eigenvalues of 
the coefficient matrix in the respective direction. Since we do not expect shocks we only consider 
a linear fourth difference artificial viscosity and not a nonlinear second difference [6], see the result 
section for more details. 

We next consider the same scheme for the preconditioned compressible inviscid equations, under 
the assumption that the entropy, S, is constant so that p = p(p)- It easier to analyze the convergence 
for the nonsymmetric form since the pressure, p, convergences and not see (10). For the 
preconditioned continuity equation we have 

02 f i 

Pt+ 0? l PX + ^ + PCl2 ( Ux + U V)J = h3 [( K \PxXx)x + {K 2 Pyyy)y\ 

Since p x = a 2 p x ^p y = a 2 p y we can rewrite the system as 


Pt + [{p u )x H" {P v )y)] ^ [{K\Pxxx)x + i^2Pyyy)y] 

"j^Pt “t" H" UUx 4* VU y H — — h V>xxx')x “h (^2 ^yyy)y] (12) 

jjPt + V t + UV X + VVy + ^ = h 3 [{K\V XXX ) X + {K 2 Vyyy)y\ 

Comparing (11) with (12) it is obvious that \i p — ► 1 as M — >0 then (11) converges to (12). 
It is crucial for both the time derivative and the artificial viscosity in the compressible continuity 
equation to be pressure based rather than density based. The preconditioning of the momentum 
equations is not important for this convergence. 

For the incompressible equations in conservative form we multiply the first equation in (11) by 
u and add it to the second and third equations. However, we do not change the artificial viscosity. 
Then 

Pt + + Vy) = h 3 [{K X p xxx ) x + ( K 2 Pyyy)y ] 

2 u 

-jpPt + U t + (u 2 + p) x + (uv)y = h 3 [{K l U xxx ) x + (K 2 Uyyy)y\ (13) 

2 'o 

-jpPt + V t + (uv) x + (v 2 + p)y = h 3 [(KlV xxx )x + (K 2 Vyyy) y ) 

For the compressible equations in conservative form we have two choices. One choice is to 
multiply the first equation in (11) by u and the second by p and add the two. The spatial derivatives 
are then in conservation form. However, the time derivative is of the form pu t rather than (pu) t 
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and the artificial viscosity terms are not in conservation form. Hence, we instead choose to apply 
the preconditioning directly to the conservative form. The resultant preconditioned compressible 
equations in conservative form is 


Pt A P 2 [{pu) x (pv) y )] — h [{KlPxxx^x + {K'2Pyyy)y] 

2 U 

+ (P U )t + {P U ^ + P)l + {P UV )y = h 3 [(Kl(pU) XXX ) X + (Ii2(pU)yyy)y] (l^) 

-jfiPt + ( pV) t + {PUV) X + (pV 2 + P)y = h 3 [(KI(PV) XXX ) X + (K 2 {pV)yyy)y) 

Note that (14) is not equivalent to (12). 

In this case we again see that (14) converges formally to (13) as M — ► 0 and p — ► 1. This is 
because the pressure is used for the time derivative and the artificial viscosity in the continuity 
equation. 

This all applies to the isentropic equations. The compressible equations for nonisentropic flow 
is more complicated and in fact there does not exist any proof of the convergence of the solution 
of the compressible equations to the solution of the incompressible equations for this case. 


6 Computational Results 


We now present a calculation for two dimensional flow around an airfoil to demonstrate the previous 
theory. As described above the discretization is based on the multistage time method coupled with 
a central difference approximation as described in [6]. 

We solve the equation in conservation form based on a hybrid set of variables of those previously 
considered. 


W t + P{F X + G y ) = AD = (K X Q 

xxx)x + (#2 Qyyy)y 


d) 





( pu ) 


* pv ^ 


p' \ 

F = 

pu 2 + p' 

, G = 

puv 

, Q = 

pu 


puv 


pv z + p 


pv 


K pH'u J 


pH'v ^ 


H' J 


where p' = p — p^, E' = c p p(T — T oo) — (p — Poo ) + P ^ U 2 ' V ^ an ^ P^' — E 1 + p 1 . We subtract the 
constants to keep the quantities in scale, see (10). 
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We choose 


B 2 

Bz 


B 4 


B iu + 


au 

w 


Biv + 


av 

w 


B\H + 


a(u 2 + v 2 ) 

w 


(3 2 = max(u 2 + + t£>)), a = l 

These equations are given for for the nondimensionalized variables. The nondimensionalization 
affects the convergence, In some codes, p and p are fixed in the far field. This implies that the 
speed of sound, a, is also bounded. As the Mach number goes to zero the pressure remains of order 
1 while the velocities go to zero. Alternatively, one can nondimensionalize so that the velocities 
are of order 1 in the far field and then the pressure and speed of sound go to infinity, unless one 
subtracts an appropriate constant, 

A typical step of a Runge-Kutta approximation is 


w (k) _ ^(o) _ ajfeAf 


D X F i*" 1 ) + DyG^ - AD ] , 


where D x and D v are spatial differencing operators, and AD represents the artificial dissipation 
terms. The dissipation terms are a blending of second and fourth differences. That is, 

AD = (Dl + D\ - D* s - D*) Q, (15) 


WUC1C . . V «. 

D lQ = A *J Qu* 

D*Q = V * [( A «+^i € !+Vi) A * V * A *] 

and A r , V* are the standard forward and backward difference operators respectively associated 
with the x direction. The variable scaling factor A is chosen as 

A *+2ij = 2 [^ Ar ^'J ( A *)*'+i.i] 

where A* and A y are proportional to the largest eigenvalues of the matrices A and B. For generalized 
coordinates x and y are replaced by £, 77 respectively. This spectral radius is now a function of the 
preconditioning. Hence, 

A* = p(PA) A v = p(PB) 

The coefficients and are adapted to the flow and are defined as follows: 




»ij = 


Pi+hi ~ 2 P»,j + Pi -1 J 


Pi+1,3 + 2 Pt,j + Pi- 1 , j 
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where p is the pressure, and the quantities and are constants to be specified. The operators 
in (15) for the y direction are defined in a similar manner. 

The second-difference dissipation term is nonlinear. Its purpose is to introduce an entropy-like 
condition and to suppress oscillations in the neighborhood of shocks. This term is small in the 
smooth portion of the flow field. The fourth- difference dissipation term is basically linear and is 
included to damp high-frequency modes and allow the scheme to approach a steady state. Only this 
term affects the linear stability of the scheme. Near shocks it is reduced to zero. For incompressible 
flow shocks can only appear in the, nonphysical, transient and so the second-difference dissipation 
is not important. To reemphasize, the preconditioning matrix multiplies the the flux terms but not 
the artificial viscosity terms. The scaling in the artificial viscosity depends on the spectral radius of 
the preconditioned matrices. If one were to use a matrix valued viscosity, [14] ,it would be related 
to the absolute value of the preconditioned Jacobian matrices. 

The boundary conditions at the far field boundary, for subsonic flow, are based on the one 
dimensional theory of characteristics in the direction normal to the boundary. The preconditioning 
changes the form of these characteristic variables. They are now given by 




where u is the component of the velocity normal to the boundary. This formulas simplify slightly 
if a = 1 and more if a = 1 + If we consider low Mach numbers then we can approximate these 

by 


Ri = u - 


Pfi' 


i?2 = « + — 

pP 


which is the same as for the incompressible case. At solid boundaries the normal momentum 
equation is used which is not affected by the preconditioning. 

The solution is advanced by the explicit Runge-Kutta method described above and without 
any residual smoothing or multigrid. We present two calculations for inviscid flow about a NACA 
0012. The first calculation is for inflow conditions M = 0.03, a = 4° . In this case we see 
that the residual asymptotes without the use of preconditioning and that the preconditioning 
dramatically increases the rate of convergence. The use of the preconditioning adds only a few 
percent to the total computational time. For viscous flows the computational time required for 
the preconditioning is negligible. In the second case we consider the same geometry but with an 
inflow of M = 0.8, a = 1.25° . In this case we also see a increased rate of convergence for the 
preconditioned case but not as dramatic as before. 

In all cases we could not allow j3 to become too small. In fact the cutoff is sufficiently large 
so that (3 is close to a constant. This has been observed in many central difference Runge-Kutta 
codes but has not been observed in the upwind code coupled with an ADI solver [3] . 


7 Conclusion 

We have considered a family of matrix preconditionings for both the incompressible and compress- 
ible fluid dynamic equations that generalize previous results. In both cases the wave speeds are 
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more equalized than for the original set of equations and so the condition number of the system is 
reduced. For the compressible equations the condition is equal to 1 at a Mach number of zero and 
increases as the Mach number increases. At M = 1 the condition number is infinite but it increases 
at a slower rate than for the physical system. 

In addition to the question of the convergence rate to a steady state we have considered the 
question of the accuracy of the numerical scheme for low Mach numbers. One can prove that for the 
partial differential equation that the compressible equations approach the incompressible equations 
as the Mach number goes to zero. For the numerical scheme this is no longer generally true and so 
the accuracy of the numerical scheme to the compressible equations decreases as the Mach number 
goes to zero. One way to improve the situation is to include the preconditioning in the Riemann 
solver, or equivalently, to account for the preconditioning in the artificial viscosity. For example, 
for low Mach numbers the scalar artificial viscosity for the continuity equation should be based on 
the pressure rather than the density. 
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A Appendix 


Let W denote the conservative variables (/?, m, n, E )*, with m = pu y n = pv , let w denote the 
primitive variables (p, and let w denote (p,u,v,Ty. Then 
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figure 1: Residual for the momentum equation. Flow about a NACA0012, 
Moo = 0.03, a = 4° with a 98x22 C grid. Graph (1) is the preconditioned 
solution and (2) is without preconditioning. 
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figure 2: Same as figure 1 but with Af* - 0.80, a 1.25 
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